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ABSTRACT 

On the basis of data from an energy- conserving 3D general relativistic MHD 
simulation, we predict the statistical character of variability in the coronal lu- 
minosity from accreting black holes. When the inner boundary of the corona is 
defined to be the electron scattering photosphere, its location depends only on 
the mass accretion rate in Eddington units m. Nearly independent of viewing an- 
gle and m, the power spectrum over the range of frequencies from approximately 
the orbital frequency at the ISCO to ~ 100 times lower is well approximated by 
a power-law with index -2, crudely consistent with the observed power spectra 
of hard X-ray fluctuations in AGN and the hard states of Galactic binary black 
holes. The underlying physical driver for variability in the light curve is variations 
in the accretion rate caused by the chaotic character of MHD turbulence, but the 
power spectrum of the coronal light output is signiflcantly steeper. Part of this 
contrast is due to the fact that the mass accretion rate can be signiflcantly mod- 
ulated by radial epicyclic motions that do not result in dissipation, and therefore 
do not drive luminosity fluctuations. The other part of this contrast is due to 
the inward decrease of the characteristic inflow time, which leads to decreasing 
radial coherence length with increasing fluctuation frequency. 

Subject headings: accretion, accretion disks — galaxies: nuclei — X-rays: binaries 
— black hole physics — MHD — radiative transfer 



1. Introduction 



In a rough manner of speaking, the light from all accreting black holes, whether those of 
stellar mass (Galactic Black Holes, or GBHs) or those residing in galactic nuclei with masses 
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10^-10^ times larger (Active Galactic Nuclei, or AGN), can be divided into a thermal and a 
coronal portion. The former corresponds to the part of the continuum spectrum with a clear 
characteristic energy (typically ~ 1 keV in GBHs, ~ 10 eV in AGN) and is thought to be 
the result of nearly- LTE emission from the surface of the accretion disk feeding the central 
black hole. The latter corresponds to the part of the continuum that extends approximately 
as a power-law from energies of order the thermal peak all the way up to ~ 100 keV and is 
thought to arise from inverse Compton scattering of seed photons provided by the thermal 
component. Both components vary. Compare d at the same (temporal) frequen cy, the coronal 
part generally varies with greater amplitude (IRemillard fc McClintockl 120061) at frequencie s 



within a few orders of magnitude of inner disk dynamical frequencies lEdelson et al.l ( 120001 ) ; 
Breedt et al.l (120091 ). It is the object of this paper to link, for the first time, the character of 
this coronal variability to heating processes directly driven by accretion dynamics. 

Coronal variability from black holes has been the subject of empirical study for several 
decades. The simplest way to characterize this phenomenon is in terms of its Fourier power 
density spectrum (PDS), -P(z^). Over several orders of magnitude in frequency, most of the 
observed fluctuation power is contained in a continuum that varies smoothly with frequency. 
It is convenient to describe this continuum in terms of its logarithmic derivative with fre- 
quency a. With the sign convention that P[y) oc z/", these slopes are commonly in the range 
-2 < a< 0. 

GBHs move among a repertory of spectral states in which the balance between ther- 
mal and coronal luminosity changes. The detailed shape of P[y) for the coronal luminos- 
ity appears to be correlated w ith the specific spectral state of a GBH (as summarized in 
Remillard fc McClintockl (120061 )). In the "low/hard" state, a tends to decrease slowly with 
increasing frequency: for example, in GROJ1655-40, a ~ for u < 1 Hz, and then gradually 
falls to less than —1 at higher frequencies, while for Cyg X-1, it decreases fr om ^ —1.7 to 
~ —2.4 over approximately the same frequency range (IRevnivtsev et al.ll2000l ). By contrast, 
in the "steep power-law" state of this system, a ^ — 1 over almost the entire observed fre- 
quency range, from 0.01 to 100 Hz. In the "thermal" state, in which the coronal luminosity 
is so weak that it dominates the spectrum only at comparatively high energies, a ~ — 1 
at low frequencies like the steep power-law state, but steepens to ~ —2 at high frequen- 
cies, somewhat resembling the low/hard state Cyg X-1, but with much smaller amplitude 
variations. Despite the coronal component's comparative weakness in the thermal state, its 
varia tions dominate those o f the thermal component, at least for photons of more than a few 
keV (jChurazov et al.ll200ll ). GBHs can also exhibit quasi-periodic oscillations (QPOs): in 
the low/hard state at ~ 1 Hz, and sometimes in the transition to the steep power-law state 
at frequencies of several hundred Hz. 
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In AGN, whose spectral state does not appear to change in the same manner, the 
power spectra may be more simply described: typically —2 <; a <; —1 below some charac- 
teristic frequency, b ut falling to <^ —2 at higher frequencies (IMarkowitz et al.l |2003| . 120071 : 
Arevalo et al.ll2008l ). QPOs, which, particularly in their lower frequency variety, are easy to 
see in the low/ hard state power spectra of GBHs, are either entirely absent or at least dif- 
ficult to dete ct (IVaughan fc Uttleyll2005l ): there is at present only one reasonable candidate 
in an AGN (IGierlihski et al.ll2008l ). The variability of AGN and GBHs can be linked (ap- 
proximately) through a simple scaling: the frequency of the roll-over in the power spectrum 
appears to be inversely proportional to the m ass, as the most naiv e theory might predict, 
albeit with a correction for the accretion rate (IMcHardy et al.ll2006l ). 



In order to understand the nearly featureless character of X-ray power spectra, many 
disparate theories have been proposed. The physical content of these theories has gradually 
increased over the years. The very first ideas were purely formal: the aperiodic natur e of the 
l ight c urves of Cyg X-1 implied to some that they were due to uncorrelated shot noise (ITerrell 
(I1972I )). Phenomenological models followed, as others supposed that power spectra that were 
power-laws in frequency resulted from a confluence of power-laws — in radius — that might 



relate the local emissivity or propensity to hot spot s to the orbital frequency (IKrolik et al. 



I991I : lAbramowicz et al.l Il99ll : iPechacek et al.l |2008| ) . These models ar e consistent with the 
observation that mar iy AGN are phase inc oherent (IKrolik et al.l Il993l ) and do not exhibit 
limit cycle behavior ( ICzerny fc Lehtol 119971 ). The fact that there can be substantial power 
at frequencies far below the orbital frequency of the inner disk coupled with the assumption 
that X-rays are predominantly emitted at small radii suggested to some that dynamics at 
large radii, where the dynamical timescales are longer , control th e low frequency behavior 
by modulating the mass accretion rate. For example, iLyubarskiil (119971 ) sought to explain 
the shape of the power spectrum in terms of fluctuations at large radius in t he disk's ratio of 



stress to pressure, the Shakura-Sunyaev a parameter. Elaborating this idea, IChurazov et al. 



(I2OOII ) suggested that such fluctuations might explain the con nection between PDS form and 
spectral state in GBHs. Alternatively, lAxelsson et al.l (120061 ) proposed that the changes in 
power spectrum might be caused by disk precession. 

None of this early work had any direct connection to physical mechanisms. Most 
importantly, in the past fifteen years we have come to realize that accretion is driven 
by MHD turbulence, an d the turbulence is stirred by the magneto-rotational instability 
(IBalbus fc Hawleyll 19981 ). To move from phenomenology to physics, models must make con- 
tact with the underlying physics of accretion. For example, studies of the long-ter m behavior 
of MH D tu rbulence might deterrn ine whether the fluctuations in a suggested by iLyubarskii 
(I1997I ) and IChurazov et al.l (I2OOII ) actually occur. It would also be highly desirable to link 
more directly fluctuations in the accretion rate with fluctuations in light output. Although 
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in the long-run the energy available for radiation is governed by the accretion rate, time- 
dependent disks may have an accretion rate that is not exactly the same at all radii, the 
heating rate may not exactly follow even the local accretion rate, time is required to gen- 
erate photons once gas is heated, and t he photons , once emitted, can take a finite time to 
make their way out of the disk. Indeed, iReig et al.l (120061 ) argue on the basis of the distinct 
variability properties of the coronal and thermal components in the soft state, and the cor- 
relation between luminosity and coronal power-law slope, that fluctuations in the accretion 
rate cannot on their own explain the observed variability in GBHs. 

More recent work has begun to follow this path, but typically has us ed proxies for the 
radiat ion rate rather than a measure of the time-varying luminosity itself. iHawley fc Krolik 
( I2OOII ) took the first step. They computed the power spectra of both the mass accretion rate 
in the plunging region and the volume-integrated magnetic stress, with the thought that one 
or the other would be a reasonable predictor of the time-dependence of the light output. The 
two power spec tra were similar, but not ide ntical, both crudely describable as power-laws 
with a ~ —1.5. lArmitage fc Reynolds! (120031 ) elaborated this approach. Assuming that the 
local emissivity follows directly from the local mass accretion rate (and not separating the 
coronal part from the thermal part), they used the vertically-integrated and azimuthally- 
averaged magnetic shear stress from 3D pseudo-Newtonian MHD simulations as a proxy for 
the local accretion rate. Placing the resulting emissivity in the disk's equatorial plane and 
assuming further that the fluid followed circular orbits, they calculated the light curves seen 
by distant observers, allowing for general relativistic ray-paths and Doppler-shifting. Even 
though the PDSs of individual radial annuli were well described by broken power-laws (ai = 
— 1 and a2 = —3.5) whose breaks were near the local orbital frequency, the superposition of 
these PDSs — because of t he radial dependence of t he em issivity — led to a = — 2 power-law 
PDSs in the total output. iMachida fc Matsumotd (120041 ) came to similar conclusions based 
on a Fourier analysis of the mass accretion rate in the plunging reg i on. M oving slightly closer 
to incorporating radiation mechanism physics, ISchnittman et al.l (120061 ) used data from 3D 
MHD simulations in full general relativity to predict model light curves and p ower spectra 
from the thermal component alone. Most recently, [Reynolds &: Millerl (120091 ) studied the 
fluctuations in a variety of dynamical quantities monitored in a 3-d pseudo-Newtonian MHD 
simulation, hoping to find an origin for QPO behavior. 

In this paper, we seek to connect dynamical calculations still more tightly to radiation. 
The tool we bring to bear on this prob l em is a new fully general relativistic 3D MHD 
simulation code (described in lNoble et al.l (120091 )). Because this code intrinsically conserves 
energy, it can self-consistently relate dynamics to heating. However, because inclusion of 
simultaneous radiation transfer is not yet feasible, we cannot provide a complete account 
of the radiation output. In particular, photon diffusion times within the disk body are so 
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long (shearing box calculations that do include rad iation transfer have shown that they are 
generally ~ 10 orbital periods: Hirose et al. ( 20061 )) that diffusion delays can substantially 
affect the time-dependence of the emerging light. Consequently, in this paper we focus on 
the variability of the luminosity from the coronal region , where optical depths are likely no 
more than order unity (see, e.g., Ilbragimov et al.l (120051 )). 



The specific example we treat is one in which the black hole rotates with spin parameter 
a/M = 0.9. Radiation is created with an emissivity in the fluid frame which depends 
on the local temperature in a way designed to give the disk a desired aspect ratio. Any 
dissipated heat is radiated on an orbital timescale; in this fashion, we attempt to make the 
time-dependence of the light output follow closely the time-dependence of heat-generation 
by such mechanisms as shocks or magnetic reconnection. The time and energy at which 
photons arrive at infinity are computed on the basis of fully general relativistic ray-tracing 
including an allowance for all travel time effects. 

The remainder of this paper is organized as follows: In Section [2] we remind the reader 
of the salient characteristics of our simulation and detail the new features of our ray-tracing 
method. In that section we also define our time-series analysis methodology. Section [3] 
presents our results, which we discuss in Section HI 



2. Methodology 



For more than a decade, MHD turbulence driven by the magnetorotation al instability 
has been recognized as the prime driver of accretion (IBalbus fc Hawleyl Il998l ). Numerical 
simulations are the most powerful tool we have for studying turbulence, and in recent years 
methods have been developed that permit simulatio ns of accretion disks over significant radial 



ranges in full 3D using g eneral relativistic dynamics (jPe Villiers &: Hawlevll2003l: IShafee et al. 



2008 



Noble et al.l 120091) . From these and analogous 2D simulations (iGammie et al 



2004 



Fragile &: Meierll2009l ). a consistent picture has emerged, despite a wide range of numerical 
algorithms and gridding schemes: Most of the accreted material flows through a dense disk 
that orbits the black hole at very nearly the angular frequency of circular orbits in the equa- 
torial plane. Within this dense disk, relatively small velocity fluctuations are superposed on 
the bulk's orbital motion. Higher in latitude, the disk becomes less dense, more magnetized, 
and more organized in both magnetic field and velocity. 



The simulation code we used to create the data discussed here was described in Noble et al. 



(120091 ). It is an intrinsically conservative ideal GRMHD code called HARM3D that accurately 
captures any gridscale numerical dissipation as heat. Numerical dissipation in many ways 
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emulates natural dissipation; when shocks collide and magnetic fields reconnect, entropy is 
created and the gas is heated. Left unchecked, the continual dissipation would make much 
of the disk unbound and lead to a progressively growing disk thickness. Both to permit 
creation of a (statistical) steady-state and to track the rate at which energy is dissipated, 
we inserted into the stress-energy conservation equation an artificial cooling function; i.e., 
this equation was given the form V^T^ = —Cu^, where V denotes a covariant derivative, 
Tj(f is the complete stress-energy tensor, C is the radiative emissivity in the fiuid frame, and 
is the fiuid four-velocity. The cooling function C = VLKtf{T /T^)^ where Vt^ is the local 
Keplerian frequency, e is the proper thermal energy density, and f{T/T^) is a continuous 
function that is zero for T/X,, < 1 and increases at higher temperatures. The local target 
temperature T* is a function of radius chosen to regulate the disk to a nearly constant aspect 
ratio H/r; in the simulation discussed here, H{r)/r ~ 0.05 — 0.12. Only gravitationally 
bound material is cooled, and (as suggested by the form of our stress-energy equation), the 
radiation is assumed to be isotropic in the fiuid frame. This simple radiation model was used 
because we are primarily interested in the bolometric emission from the disk and wish to 
apply it to a wide variety of black hole sy stems. A more model-dependent cooling function 



could also be used (IFragile fc Meierl 120091 ) . but it would be computationally more expensive 



and would also require choosing both a specific black hole mass and an accretion rate. 

Our numerical domain was divided into 192 x 192 x 64 cells in the radial, poloidal, and 
azimuthal directions respectively, with r G [1.28, 120]M, 9 e [O.OStt, 0.957r], G [O,7r/20. 
The radial discretization is logarithmic — Ar oc r — to resolve finer features at smaller radius. 
The azimuthal resolution is constant, and the poloidal discretization is rarefied at the poles 
and concentrated at the equator. 

The pressure maximum of the initial distribution — at r = 25M — sets the location within 
which a well defined accretion fiow exists. The disk reaches an infiow steady-state for r < 
14M over the period t = [7000M, 15000M]; we examine only this epoch here. For reference, 
the orbital period at radius r is Torhir) = 3.1 x 10"^ (M/IOMq) {r/Mf^"^ + a/M s. The 
span At = 8000M represents approximately 287 orbital periods at the innermost stable 
circular orbit (ISCO) and 10 orbital periods at the initial pressure maximum. For our black 
hole spin parameter, risco = 2.32M and the horizon is located at rhor = 1.44M. The disk's 
rest-mass density p, 4-velocity and cooling function C evaluated at all grid points are 
written to disk every 20M in time. We use this data as input to our radiation transfer 
procedure to create light curves. Any emission outside r = 25M is ignored. 



^Note that throughout this paper we use geometrized units with G — c = 1 unless mentioned otherwise; 
distances and times are therefore scaled to the mass of the black hole M. 
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Because the focus of this paper is time-variabihty properties, we point out that our 
method has two hmitations that affect the shortest timescalcs. Samphng at intervals of 20M 
means that no frequencies higher than 1/(40M) can be probed; our Nyquist frequency is 0.7 
times the orbital frequency at the ISCO. The other limitation comes from our cooling rate. 
Because the characteristic cooling rate is ~ Qk, heating fluctuations on timescales shorter 
than ~ flj} cannot be translated into equally rapid emission fluctuations, even though 
some cooling mechanisms, notably inverse Compton scattering, can often have cooling rates 
considerably faster than ~ D,k- In sum, we cannot present results on frequencies above 
~ O.lvisco, and the form of our cooling function potentially suppresses some fluctuation 
power at the higher frequency end of the range we do discuss. 



2.1. Radiation Transfer 

Within the simulation, we do not consider any interaction between the emitted radiation 
and the material. However, more realistically, there is always some opacity, and in most 
circumstances the dominant opacity in the material near a black hole is electron scattering. 
This opacity leads to a natural division of the radiation in two parts: that emitted inside or 
outside the photosphere. Within the photosphere, scattering can add substantially to the 
time a photon can take to reach the outside, washing out fluctuations in intrinsic emissivity; 
outside the photosphere, of course, scattering has very little effect on photon escape time. In 
addition, photons deriving their energy from dissipation inside and outside the photosphere 
can be distinguished spectrally: Inside the photosphere, thermalization is strong, and the 
local spectrum should be approximately black body, at a temperature ~ 10 eV in AGN, 
~ 1 keV in GBHs. By contrast, outside the photosphere, much lower gas densities and much 
higher ratios of heating density to mass density lead to much higher temperatures, and the 
primary emission mechanism is inverse Compton scattering, so that the radiated spectrum 
is characteristically a power-law extending well into the hard X-ray regime. In order to 
make a realistic estimate of the light curve directly from the simulation's emissivity data, 
we therefore restrict our efforts to the coronal hard X-ray emission, whose source is near or 
outside the scattering photosphere. 

To locate that photosphere requires a calculation of the opacity, yet its magnitude is 
not deflned in code-units because the simulation requires no absolute density scale. Instead, 
we determine it after the fact by the following procedure: We distinguish quantities in code- 
units from quantities in physical units by attaching a subscript c to the former, and leaving 
the latter unlabeled. If a fraction 77 of the rest-mass of accretion were transformed into 
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luminosity at infinity, it would be 

L^vJ ded(l)V^p,u'{p/p,){GM/cYc^ = r)Mc{p/pc){GM/c^fc^ (1) 

because the unit of lengtli is GM/c^ if G = c = 1, and u^, wlien measured in units of c, 
is dimensionless. Here, g is tlie determinant of tlie metric. Normalizing the luminosity to 
the Eddington luminosity Le, we find that the relation between physical density and code 
density is 

p p^ = : — , 2 

' ktGMMc vLe ^ ^ 

where kt is the electron scattering opacity per unit mass, and — 0.0177 is the time- 
averaged rest-mass accretion rate in code units. By fortunate coincidence, optical depths 
depend only on L/{'r]LE), which we abbreviate as m, because the unit of length is oc M. 

Because our accretion flow is far from spherically symmetric, the location of the photo- 
sphere is a function of the observer's position. We imagine, then, that numerous "cameras" 
are placed on a grid in polar angle ?9 and azimuthal angle (/? on a very large sphere (radius 
lO^M) centered on the black hole. From each camera, we define a bundle of geodesies that 
run through the problem volume. These are parameterized by an affine parameter A normal- 
ized so that an observer in the local fluid frame would measure the differential length along 
a ray as 

ds = udX, (3) 

where u is the frequency of the photon as measured by that observer. If N'^ = dx^ / dX is the 
4- vector tangent to the null ray then 

1/ = ^ (4) 

z 

where z is the redshift factor between the local fluid frame observer and the camera frame: 

In the numerator of this ratio, the 4- velocity is that of the camera; in the denominator, it is 
that of the fluid at some point along the ray. We then integrate the optical depth 

dr = pKTi^dX (6) 

along these geodesies in order to determine the location of the photospheric surface for that 
camera. The photosphere surface is defined to lie at a constant r = Tq, which we set to unity, 
i.e.. To = 1. 
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Once the location of the photosphere is determined, we integrate the emissivity along 
these geodesies from the photosphere out to the camera; we assume no scattering takes place 
along these rays: 



where 1^ is the specific intensity and is the fluid-frame emissivity, given by: 

Jj. = ^ {ZV - i^cam) ■ (8) 

Integrating over all i^cam to find the bolometric luminosity is equivalent to setting v^axa — 1 
and V — 1/z] the latter procedure is done in practice. To set the units of the observed 
luminosity, we note that the units of power density are the units of energy density (pc^) 
divided by the unit of time (GM/c^). The end result is 

L = ^ . 9 

KT{GMfMc 

However, these units are also unnecessary because all our results for variability will be shown 
in fractional terms, relative to the mean luminosity. 

We are therefore left with three parameters to explore: m, 'd and ip. We vary rh from a 
value low enough that the entire flow is optically thin up to the Eddington limit: 

m e {0.001,0.003,0.01,0.03,0.1,0.3,1.} . (10) 

The simulation should not be biased toward any particular pole, so we sample '& only over 
one hemisphere, uniformly in sin d: 

1? e {5°, 17°, 29°, 41°, 53°, 65°, 77°, 89°} . (11) 

Similarly, the physics of our accretion disk has no special azimuthal orientation, so any 
observed dependence of the light curves on must be only statistical fluctuations. However, 
our simulation domain spans only the flrst quadrant in azimuth, from to 7r/2. 

To cope with this limitation, we remap the density and velocity data into the other 
quadrants, but not the emissivity. By doing so, we can compute the portion of the light 
reaching inflnity from this quadrant alone with a proper allowance for optical depth effects 
in all directions. In principle, there are four different ways we might have placed the radi- 
ating quadrant with respect to the quadrants having only opacity. Prom the expectation of 
azimuthal symmetry, it then follows that a full description of the statistical character of the 
light curve can be obtained from viewing this quadrant from only four azimuthal directions, 
which we choose as 

TT 37r 1 



(^e<!0,-,7r,-| . (12) 
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This piecewise construction of the composite image is illustrated in Figures [D - El Images 
from different (y? are shown in Figure (H and their sum is plotted in Figure [2l We find that 
frame-dragging of photons can cause emission from one quadrant of the simulation domain 
to spread into all quadrants of the image plane. The total flux of a quadrant varies with 
since certain values orient the disk's orbital velocity closer to the line of sight. Figures [1] 
and [2] will be helpful references in our later discussion about how the artificial azimuthal 
symmetry condition influences our variability predictions in Section 13.41 

Tracking the light through the simulation data is complicated by the fact that space- 
time's curvature means that a set of photons reaching an observer at one instant may orig- 



i nate from the accretion flow over a wide range of coordinate times. lArmitage fc Reynolds 



(120031 ) included the effect of time delays, but used an infinitesimally thin emission region 
which intersected each geodesic only once. This meant that they did not have to propa- 
gate the light rays through the simulation data, as we must because our emission region is 
extended. 

There are many ways to go about the bookkeeping inherent to this problem, but dif- 
ferent algorithms place vastly different demands on computer memory. Ray-tracing as the 
simulation runs is expensive, memory intensive, and not amenable to exploration since ad- 
justments made to the ray-tracing scheme would require rerunning the simulatior|^. For this 
reason, we adopted a post-processing procedure, which means that our time resolution is 
set by the rate at which we output simulation data. As a ray traverses spacetime, we use 
quad-linear (linear in space and time) interpolation to determine the necessary quantities 
along the ray's path. This interpolation is done with pairs of simulation data slices at a 
time in order to reduce the calculation's memory footprint. Rays are organized into stacks 
of snapshots, each representing a batch of rays distributed over the image plane that reach 
the observer synchronously. Each snapshot then depends on a finite span of simulation data, 
or rather, a given time slice of simulation data influences a sequence of snapshots. 



In iNoble et al.l (120091 ). we spent much effort to ensure our ray tracing calculation was 
converged with respect to the number of light rays used per image; following the analogy with 
a pinhole camera, we will refer to them as pixels in our camera. We found best performance 
using a nonuniform pixelation that approximates a projection of the simulation's grid onto 
the image plane. Using the radial profile of flux at infinity dF/dr as our convergence criterion, 
we found that iVpix = 255^ nonuniform pixels and iVpix = 1200^ uniform pixels result in 
approximately the same level of accuracy (largest relative error over r is 5%) when compared 
to a calculation with Npi^ = 1500^ uniform pixels. 



^See De Villiers ( 20081 ) for an algorithm that performs the ray-tracing in situ. 
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2.2. Time Series Analysis 

The radiative transfer calculation yields a function F {i}, t, m) which represents the 
bolometric flux for an observer at constant radius (r = lO^GM/c^) and time t. In the 
following analysis unless explicitly stated, we integrate F over to yield a light curve for a 
full 271 disk. Our goal is then to understand the temporal character (e.g., power spectrum, 
variance) of this function with respect to the remaining parameters (i.e. d,rn). 

Given a function defined on a set of discrete times F{tn) (where n = 0,...,A^ — 1), 
we define its discrete Fourier transform F{i>) at frequency u as 

N-l 
n=0 

where u G {1/T, . . . , N/2T}, and T = t^^i — to- As remarked earlier, our Nyquist limit is 
t'max = (40M)~^ ~ O.Tuisco, where z/isco is the orbital frequency at the ISCO. 

We define the normalized PDS of the light curve F{t) such that it is a measure of the 
fractional variance per unit frequency: 

(14) 

where F is the time average of F{t) n Since F {i') is a complex number, we can represent it 
in terms of its magnitude and phase: F (u) = A (u) e*'^*^^^. 

Observed light curve power spectra are often described in terms of a best-fitting power- 
law; as we will see, our results resemble power-laws at about the same level of accuracy. 
Determining the best-fit slope to the power spectra of our light curves, however, is somewhat 
arbitrary because we do not know the "error distribution" of our data (the relevant ensemble 
for us would be a large set of simulations begun with slightly different initial conditions). 
Because our goal is only quahtative description, we choose a very simple approach: equal 
standard error at each frequency point. That is, to find the best-fit power-law slope for a 



Fiv] 



^ These power s pectr a acc urately descri be our finite duration simulation. However, as pointed out by 
Deeter fc Bovnton ( 1982h and Deetei (1984), power spectra with slopes steeper than —2 may be subject to 



artificial leakage of fluctuation power from low frequencies to high. This problem can be significant when the 
power spectrum is substantially steeper than —2 and extends to frequencies significantly lower than those 
probed by the experiment (i.e., much lower than the inverse of its duration). Because the slope we measure 
is only slightly more negative than —2 over most of the parameter space of interest, we believe that this 
artifact may have only limited effect, but only much longer simulations can definitively answer this question. 
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given power spectrum, we minimize the quantity 



Xeff 



J2 [Pi^i) - CK] 



(15) 



with respect to C and a. This form tends to weight most heavily the highest decade of 
frequencies because the frequency points are separated by a constant Az/ = (SOOOM)^^ 

A key physical question we wish to explore is: how well correlated is the variability 
at different radii? Answering this question will help us understand the light curve ' s po\Y er 
spectrum and provide a context for comparison to other models (e.g. lLyubarskiil (119971 )). 
Correlations are performed by first decomposing the light curve's PDS into partitions. Let 
there be N partitions, each with a different light curve F„ (t) so that our total light curve 
is their simple sum: F [t) = Yln=iPn (t). The total power spectrum P(z^) can therefore be 
expressed as a sum of the partitions' PDSs (Sa) plus a sum that depends on how well the 
different modes are correlated (S'f,): 



P(iy) = ^ F 



2T 
F2 



2T 
F2 



m 

^ ^2 + 2 ^ ^ AnA^ COS {A^nrr 



n 
=12 



m n>m 



(16) 



;i7) 



F AT 



F2 

n 

Sa + Sb- 



m n>m 



(19) 



Here Ai^nm = 4'n — 4'm is the difference in phase at frequency u between two partitions. 
Note that even though P(z/) is independent of our partition scheme, the relative sizes of 
Sa and Sb are not. If the An are all of similar magnitude, then Sa/ Sb — >■ as — > oo 
{Sb ^ N {N — 1) while Sa ~ A^) and Sb/Sa — > as A^ — >• 1 (there is only one partition and 
one signal is perfectly coherent with itself). If the partitions are perfectly incoherent, then 
P ~ Conversely, if they are perfectly coherent, then P Sb for N > 2. 



3. Results 

3.1. Light Curves and Power Spectra: Dependence on Accretion Rate and 

Inclination 

Sample light curves and their corresponding power spectra can be seen in Figures [3] and 
m the former displaying how they change with viewing angle at fixed m, the latter how 
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they change with accretion rate m at fixed Changing inchnation does relatively little to 
alter long time-scale fiuctuations, but can lead to differences on short time-scales. On the 
other hand, changing the opacity can lead to substantial differences in the light curve even 
while the viewing angle remains the same. Remarkably, however, the gross shape of the 
power spectrum is almost invariant to both sorts of changes: the best-fit power-law index is 
~ —2 for all but the highest accretion rates and inclinations (Fig. [5]). 

The strongest effects infiuencing the inclination dependence of variations are relativistic 
beaming and boosting, which become more important as the orbital velocity becomes larger 
and more nearly parallel to the outgoing geodesies. They therefore have the greatest effect 
on radiation issuing from the smallest radii when viewed at high inclination. Because those 
same inner radii have the highest dynamical frequency, one might then expect a boost in the 
high-frequency portion of the power spectrum at large i). In relative terms, this does not 
occur — as we have seen, the slope of the power spectrum depends only weakly on inclination, 
except when m is quite large (see further discussion later in this subsection). Nonetheless, 
although the relative variance changes little with inclination, the absolute variance, as well 
as the absolute luminosity, does in crease when the disk becoraes more edge-on, as has also 



been seen in previous calculations (lArmitage fc Reynolds! l2003l : ISchnittman et al.l 120061 ). 



Because we explore only relative variability, the absolute luminosity's proportionality 
to m is irrelevant to our discussion. The accretion rate influences the light curves in our 
calculations only by setting the opacity scale. The accretion rate is therefore degenerate with 
our choice of Tq, and we can speak equivalently in terms of accretion rate or optical depth. 

When the opacity is dominated by electron scattering, the disk is completely transparent 
for accretion rates m = 0.001 or lower. Increasing m moves the photosphere farther from 
the disk's midplane, and emission from high latitudes becomes more dominant because our 
disk follows a nearly constant H/r profile. At the same time, increasing rh leads to a relative 
suppression of light from outer radii because the disk surface density, and hence its optical 
depth, increases rapidly outward. For this reason, the largest accretion rates select out 
fiuctuations from the innermost and uppermost regions of the (bound) accretion fiow. 

This pruning of the coronal volume with increasing rh is the most likely explanation for 
the fact that the relative variance of the light curves monotonically increases with accretion 
rate, from 0.04 at m = 0.001 to 0.09 at m = 1. As the region above the photosphere shrinks 
in radial and vertical extent with m, it contains fewer independently-fluctuating volumes, so 
that their summed emission has larger fractional fluctuations. 

Increasing m also leads to greater obscuration of high inclination observers' views of 
the inner disk. It is this effect that explains the steepening of the best-fit power-law in the 
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high rh and i} corner of Figure [51 Since we restrict the emission to r < 25M while the disk 
matter extends out to r = 120M, sufficiently large accretion rate and inclination angle can 
lead to complete obscuration of the emission region. The radius within which all emission is 
obscured, Ro, is consequently an increasing function of inclination and accretion rate. Curves 
showing how Ro depends on rh and are also shown in Figure [51 Only when > 60° is this 
obscuration effect significant. Because the location of the obscured region is sensitive to the 
geometry of the disk — which is artificially tuned to have a near constant aspect ratio — the 
steepening of the fluctuation power spectrum due to obscuration may be artificial. 



3.2. The Optically Thin Limit: Origin of the Variability 

The most obvious explanation for variability in radiative output is variability in the ac- 
cretion rate. Let us ffist examine the relationship between the local emissivity and accretion 
rate to see if this is indeed true for our simulation. In Figure [SI we compare the accretion rate 
and emissivity at r = 3.5M to the light curve (integrated over the entire simulation) mea- 
sured face-on (i) = 5°). We choose to compare behavior at r = 3.5M to the total light curve 
because it is the radius of the brightest annulus, and should therefore be a major contribu- 
tor to the composite light curve. A nearly-polar viewing angle minimizes relativistic effects, 
simplifying the comparison of the observed light curve to local emissivity. As expected, the 
disk-integrated light curve follows the same large amplitude, long timescale fluctuations seen 
in the accretion rate and emissivity at r = 3.5M. However, it lacks the short timescale 
variability of the local emissivity and accretion rate. The same effect appears, of course, in 
the power spectra. At this radius (and at most of the others in the steady-state portion of 
the disk), the accretion rate and emissivity power spectra are approximate power-laws with 
exponents ~ —1 and ~ —1.5, respectively, significantly shallower than the total light curve 
power spectrum, for which the overall slope is —2. 

In order to elucidate why fluctuations in the local properties have more high frequency 
power than the total light curve, and to understand better to what degree the accretion 
rate drives the radiation, we examine the relationships between the power spectra of M(r, t), 
C{r,t), the flux from r as it is observed on the polar axis at infinity (dF/dr), and the total 
flux F{t) in the polar direction. We have already seen that |Mp and |£p are similar but 
not identical. A closer comparison of these two power spectra may be seen in Figure U\ 
which shows the ratio of the emissivity's power spectrum to that of the accretion rate as a 
function of radius and frequency. In much of the diagram, the ratio is near unity, but there 
is a depression in the ratio along a track whose frequency falls with increasing radius. This 
dip seems to be due to an excess in fluctuation power in the accretion rate; in the case of 
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r = 3.5M, a bump in |Mp can be clearly seen at z/ ^ (0.2-0.4) z/jsco (Fig. E]). 

The origin of this excess can be identified by comparing its track in radius-frequency 
space with the radial dependence of several significant dynamical frequencies: the orbital 
frequency, the vertical epicyclic frequency, and the radial epicyclic frequency. As can be 
clearly seen in Figure [3, the radial epicyclic frequency — but not the others — follows closely 
the centerline of the dip, suggesting that radial epicyclic modes modulate the accretion rate 
without giving rise to emission — i.e. the oscillations are either not dissipated significantly, 
not dissipated locally, or both. A si milar feature in th e powe r spectrum of radial velocity 



as a function of radius was noted by [Reynolds fc Millerl (120091 ) in their data from a pseudo 



Newtonian global disk simulation. The accretion rate also has more fiuctuation power than 
the emissivity for r > risco and v ~ 0.3z/igco- We do not understand the origin of this excess. 
The net result of both excesses, however, is to make |Mp a fiatter function of frequency than 
|£p at most radii and also to create deviations from power-law behavior in the local accretion 
rate power spectra. 

We next study how closely fiuctuations in Or) are mirrored in dF/dr. To do so, we 
look at the ratio 



n 



J C{r,t)dt 



f^dt 



(20) 



This is the ratio of the two normalized power spectra as a function of radius and frequency 
(Fig. [8]). We find that TZ is evenly distributed about unity, with deviations that rarely exceed 
a factor of 2 in either direction. Thus, the emissivity at r predicts dF{r)/dr a.t 'd = quite 
well. 

We focus next on how the individual annular contributions to the fiux dF/ dr sum to the 
total fiux F{t). One clue is given by the fact that the fractional variances of C{r = 3.5M) and 
M(r = 3.5M) are rather similar, 0.152 and 0.175, respectively, while the fractional variance 
of F is rather smaller, 0.029. We now understand that the emissivity follows the variability 
of the accretion rate (but with certain exceptions like those associated with radial epicyclic 
motions) and dF/dr varies like the emissivity. Why, though, does the total fiux have such 
a small relative variance, and how can a set of oscillators (disk annuli) with power spectra 
that are a ~ — 1.5 power-laws integrate to have a composite PDS with a ~ —2? 

In the language of Section 12.21 the annuli can be thought of as partitions with their own 
individual light curves. Since there are a large number of annuli or partitions, Sh > Sa unless 
there is a dramatically low degree of phase- coherence between the different radial segments. 
If all the annuli were perfectly coherent, Aipnm = Vn, m, P c::^ Sb and the light curve would 
have a a ~ —1.5 power-law power spectrum with a larger variance. On the other hand, if 
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all the annuli were completely incoherent, the total flux power spectrum would still be a 
power-law with a ~ —1.5, but with a smaller variance. The only way to steepen the slope 
of the spectrum is for the degree of coherence to decline with increasing frequency. 

The level of coherence in the variability of dF/dr is illustrated in Figure [9l where 
we plot ip{v,r) over the lower half of our frequency range. At almost all radii, ip^u^r) is 
incoherent in frequency (negligible correlation lengths in z/), but at fixed frequency, there 
can be significant coherence in radius. The phases are sufficiently coherent between different 
annuli that ^ Sa-, but their correlations follow no simple pattern. Different frequencies 
show different radial coherence patterns, making it impossible to state that radius r varies 
coherently with radius r'; rather, one can only say that certain modes at r are coherent with 
those at r'. 

The white dashed curve in Figure [3 shows the inflow rate Vi^^iow as a function of radius, 
which we have defined to be the mass-weighted mean radial velocity of bound material 
divided by the local radial coordinate: 

f, , dO dd) dt\/^ pu^ , , 

Jbound ded(Pdt^p 

A fiuid element is considered bound if hut > — 1, and h is the fiuid element's specific enthalpy. 
The time integral is performed over our standard epoch of t = [7000M, 15000M]. For 7M ^ 
r <^ 20M, we find that I'mHowir) ~ [28Torh{r)]^^ . At smaller radii, the infiow accelerates 
until near the ISCO and in the plunging region UmHowir) ~ Qk- Regions to the left of 
this curve are clearly more coherent than those to the right. That this should be so is not 
too surpris i ng, gi ven the ultimate dependence of energy release on mass infiow. Indeed, 



Lyubarskiil (119971 ) proposed that the inner disk's low frequency variability can be entirely 
explained by variations spawned at larger radii (by fiuctuations in the stress to pressure ratio) 
that are then advected inward with the accretion fiow. What is demonstrated in this phase 
picture is that fiuctuations lower than the local infiow rate do indeed propagate coherently 
inward, whatever their initial source. However, over much of the range of frequencies studied 
here, this criterion can be satisfied only near the ISCO and in the plunging region itself. At 
these higher frequencies (which, as we shall see in the next section, are often the object of 
most observational study), no such regular propagation pattern can be discerned. 

Returning to the question of why the power spectrum of the total fiux is steeper than 
that of the fiux radiated by individual annuli, we now see that this can be explained by the 
diminution of the coherent radial range with increasing frequency shown in Figure O 
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3.3. Finite Speed of Light 

Our calculation for the first time correctly accounts for time delays while ray-tracing 3D 
GRMHD simulation data. We would therefore like to quantify how our results are affected 
by inclusion of this effect. The light curves and power spectra from calculations with and 
without time delays are shown in Figure [101 The light curves are identical except that the 
time delay calculation shows slightly less short timescale variation. This fact is illustrated 
more clearly in the power spectra panel of this figure, which clearly shows that the fluctuation 
power at high frequencies is diminished when one includes time delays. 

This contrast is easily understood. Delay effects can diminish coherence in the received 
signal when a region whose light crossing time is At varies coherently on timescales shorter 
than At. On the other hand, enhancement of coherence by delay effects would require 
remarkable contrivance because spatial and temporal fluctuations in the turbulence would 
have to be correlated with the ray trajectories for particular observers. Consequently, photon 
time delays in general decrease the fluctuation power. The depression of the fluctuation power 
is confined to the highest frequencies because maintenance of emissivity coherence requires a 
coordinating signal propagating across the region, but all signals, whether conveyed in bulk 
fluid motion or by some wave mode, are limited to traveling no faster than c. It follows 
that, for light travel time effects to suppress variability, the coordinating signals must be 
relativistic. In the context of an accretion flow, relativistic signals are largely confined to 
the innermost regions, which dominate the generation of high frequency fluctuations. 

Because the time delay effect depends on the light's path through the material and the 
local velocity of the fluid, one expects it to depend on i) and m. We characterize its trend 
over parameter space in Figure [TT|, where the difference in power-law exponents between 
the calculation with time delays and that without time delays is plotted. In all cases, the 
time delay calculation yields a steeper PDS. The contrast depends most strongly on accretion 
rate, in the sense that it diminishes as the disk becomes more opaque; this trend is consistent 
with the observation that as m increases, the inner portion of the disk becomes progressively 
more obscured and contributes less to the power spectrum. Larger t? produces slightly larger 
deviations between the two methods. As the inclination angle increases, photon rays become 
more nearly parallel to the disk's orbital velocity. For those fluid elements with relativistic 
velocities, the result is that fluid elements' worldlines move closer to the lightcone, leading 
to somewhat greater coherence of emissivity along the rays. 
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3.4. Azimuthal Symmetry Condition 

In order to save computational resources, we assumed in our simulation that the disk is 
periodic in cf) on intervals equal to 7r/2. Any modes longer than this are not included in our 
calculation, and the symmetry condition artificially gives rise to correlations at this scale. 
The question then arises: are the absent modes and artificial symmetry important to our 
prediction? 

Even though we cannot rigorously evaluate the constraint's ramifications without re- 
peating our analysis for a run on the entire dom ain, some insight might come from a similar 



simulation that used the full azimuthal extent (ISchnittman et al.ll2006l ). Their calculation 
used a numerical code that ina dvertently preserves near-co nstant aspect ratio by failing to 
capture all dissipation as heat JPe Vilhers fc HawlevI 120031 ) . Even though their numerical 



techniques were different and no explicit cooling was used, our calculations share nearly 
identical initial conditions (besides the full azimuthal extent) and yield similar disk thick- 
nesses. Since the disk's thickness dictates the poloidal size of turbulent eddies in the bulk, 
we may expect that the characteristics of their correlations in will be applicable to our 
system. They found that the surface density's dominant azimuthal correlation lengthscale 
is approximately 0.47r, suggesting that our grid may be large enough to include the most 
important modes. 

To quantify the systematic effect of the 7r/2 periodicity, we can employ the partition 
formalism previously introduced. For the purpose of this discussion, it is convenient to 
label the quadrants by their panel labels in Figured], i.e., a, b, c, and d. The quadrants 
are distinguished by the sign of their mean line-of-sight velocity (receding or approaching) 
and their position (front or back). Quadrants a and c are approaching, whereas b and d 
are receding; quadrants a and b are in back of the black hole, whereas c and d are in front. 
Because relativistic effects dominate obscuration effects in determining the characters of their 
light curves, it is easiest to think of the system in the optically thin limit. If only special 
relativistic effects applied, the two receding quadrants would produce identical light curves, 
as would the two approaching quadrants. However, general relativistic frame- dragging and 
light deflection complicate the story. For example, light in a particular direction in the fluid 
frame can be wrapped around the black hole and escape at a completely different angle. 
Although the deflection angle is large only very near the black hole, most of the light is 
produced at these same radii, so it can be a very important effect. For instance, looking at 
Figure [1] we flnd that the most intense part of quadrant b is located on the opposite side of 
the black hole in the image plane because the brightest light — that which is emitted along 
the orbital velocity — has been bent around the black hole and focused toward the observer. 
This phenomenon means that quadrant b is more like an "approaching" quadrant than a 
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receding one. 

With these points in mind, we can now explain what controls the distinctions in flux 
and power spectrum between the different quadrants. For face-on views, they all contribute 
identically to the light curve; as the viewing angle moves off-axis, special relativistic beaming 
and boosting enhances the approaching sides, while general relativistic light bending and 
frame-dragging enhance the back sides. The result is that over most of rh-^ parameter 
space, quadrant a is the brightest (both approaching and in back of the black hole), d is the 
faintest (both receding and in front), and b and c are similar to one another (b is receding but 
in back; c is approaching but in front). The maximum flux contrast between the brightest 
and dimmest quadrants never exceeds a factor of 5. The slopes of their power spectra 
follow the same trend seen in flux: q„ ^ ai, c:^ > on average, with no spectral slope 
faUing outside the range —2.4 < a < —1.8. However, relative to quadrant c, the quadrant 
b becomes brighter and its PDS flatter as ■»? increases. In the summed light curve, the 
contrasting effects largely cancel one another, so that the spectral slope of the composite 
PDS can be described by a simple average of the quadrants' individual power-law exponents. 

In our calculation, the quadrants are precisely coherent at all frequencies when viewed 
exactly face-on. As the inclination angle grows, they begin to become incoherent at the 
highest frequencies, but even for 'd = n/2, the range of incoherent frequencies is still quite 
limited. The reason for this behavior is that our symmetry condition makes their emissivity 
precisely coherent, so such incoherence as exists is entirely due to time-delay effects; as just 
discussed, they arc small except at the highest frequencies. Thus, if the absolute power 
spectrum from a single quadrant (before Doppler adjustments and obscuration effects) is 
A'^^u), our total power spectrum is IGA^^u) when viewed on-axis, and when viewed off-axis 
has essentially identical power at low frequencies, but slightly less at high. 

By contrast, in a full 27r simulation we expect that the emissivitics of the quadrants 
would have very similar power spectra to the emissivity we calculate, but be completely 
incoherent if azimuthal correlations extend only over angles ~ OAtt. The same repertory of 
relativistic effects, both special and general, will still apply, but we expect that they will 
similarly cancel in sum. Thus, a total flux power spectrum ~ 4A^(i/) should result, as only 
the Sa term contributes. In other words, if this reasoning holds, the shape of the power 
spectrum observed from a full 2n disk would be quite similar to what we compute, but its 
amplitude would be lower by about a factor of 4. 
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4. Discussion and Conclusion 

In this paper, we have presented a new, more physical method for estimating the tem- 
poral variability of radiation from the optically thin ("coronal") regions of 3D GRMHD 
simulations. For the purpose of investigating variability, it is necessary at this stage to sep- 
arate optically thin regions from optically thick because present-day global disk simulation 
codes do not have the capacity to solve the transfer problem simultaneously with the dynam- 
ics, and diffusion through optically thick regions materially alters the character of variability. 
The key improvement over previous calculations is the use of data from an energy-conserving 
code with precise control of the disk's thermodynamics. In addition, we have shown how 
to include photon travel time delays, although they have a relatively small impact on the 
results shown here. We separated coronal emission from disk emission by integrating the 
fluid emissivity from the scattering photosphere outward; the location of the photosphere 
moves in a manner controlled by the nominal accretion rate in Eddington units, rfi. Because 
our density data — which determines the optical depth — was written only every 20M in time 
over a period of 8000M, we were constrained to explore the disk's power spectrum only over 
the frequency range u e [3.5 x 10~^, 0.7] i^isco- 

We found that the power spectrum of the observed flux's fluctuations in this frequency 
band is described well by a featureless power-law with index a ~ —2 for essentially all 
optical depths (or, in this formalism, accretion rates) and inclination angles. Although most 
of the fluctuation power has its physical origin in accretion rate fluctuations, the slope of 
its power spectrum is steeper by ~ 1 than the slope of the accretion rate's power spectrum. 
Two separate effects combine to create this steepening: there is high frequency power in 
the accretion rate due to radial epicyclic motions that do not contribute to variations in the 
emissivity; and the radial coherence of different frequency modes declines with increasing 
frequency. Thus, the power spectrum of accretion rate fluctuations is not a good proxy 
for fluctuations in the coronal light. Because photon diffusion damping of high-frequency 
emissivity fluctuations will likely steepen the power spectrum of the thermal luminosity, we 
expect that the same will be true for the light curve of the thermal component. 

Relativistic beaming and boosting cause the variance of the light curve to increase with 
inclination angle, but do not materially change the shape of the power spectrum. The reason 
for this perhaps counter-intuitive result is that Doppler effects flatten the power spectrum of 
the approaching segment of the disk and steepen the power spectrum of the receding segment 
so that the two changes compensate for one another. 

Time-delay effects steepen the observed power spectrum at the highest frequencies. If 
we had saved data from this simulation at intervals shorter than 20M, we expect that this 
effect would have been increasingly important at the higher frequencies that would then have 
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been accessible. For this reason, in any future work using simulation data to predict light 
curves on short timescales, we strongly encourage proper accounting for photon travel times. 

Our results span the frequency range ~ 2.5Hz — 500Hz if scaled to a M = IOMq black 
hole and ~ 2.5 x 10^'^Hz - 5 x lO^^Hz for a M = 10^ Mq black hole. Low frequency QPOs, 
which generally occur in GBHs at ~ 1 Hz, would therefore be at best only marginally 
detectable in our data. High frequency QPOs, which are sometimes seen at ~ 200-300 Hz, 
might in principle have appeared, but we see no evidence for any. On the other hand, 
they appear in real black hole systems only in association with the transition to the steep 
power law state. Although our simulation code very accurately conserves energy, the cooling 
function we employ is no more than a toy-model. A more complete description of radiative 
cooling will be necessary to understand spectral state transitions, and that might be a 
prerequisite for understanding high frequency QPOs as well. 

Our optically thin limit (m = 0.001) phenomenologically resembles the hard state of 
Galactic black hole binaries in the sense that our definition of "coronal" eliminates any 
optically-thick thermal disk in the inner part of the accretion flow. Intriguingly, the power- 
law slope that we consistently find {a ~ —2) is crudely consistent with the mean slope of 
the power spectrum measur e d in C yg X-1 in the range 1-500 Hz: steepening from —1.7 
to -2.4 



Revnivtsev et all ((20001) 



For higher accretion rates, our corona is restricted to the outer layers of the flow, more in 
keeping with what is often imagined for AGN. A power-la w slope ^ —2 is al so very roughly 



consistent with observations of these objects. For example, iMarkowita (120091 ) shows that the 



slope of the power spectrum in I C 4329A steepens froin ~ — 1 to ^ —2 across the frequency 



range 10~^-10~^ Hz. Similarly, iMarkowitz et al.l (120071 ) find that the power spectrum of 
Mrk 766 steepens from ~ —1.5 to ~ —3 from ~ 3 x 10~^ Hz to ~ 10~^ Hz. In this latter 
case, Markowitz et al. estimate that the central black hole mass may be only ~ 10^-10 ''M©, 
which would place our simulated frequency range roughly coincident with the observed band. 

As already mentioned, the decrease in radial coherence length with increasing frequency 
steepens the power spectrum of the aggregate light curve relative to the power spectrum of 
the local emissivity. In our very approximate treatment, we described the result in terms of 
a new, steeper power-law. A more careful and complete treatment might improve upon this 
description. In particular, the factor that controls the radial coherence length is whether 
the fluctuation frequency is larger or smaller than the local inflow rate. It is the outward 
decrease of the inflow rate that leads to higher power at lower frequencies by stretching the 
range of radial coherence. However, at sufficiently low frequencies, greater radial coherence 
does not add appreciably to the power spectrum because material at larger radius does not 
contribute much to the luminosity. At frequencies lower than the inflow rate at the radius 
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within which most of the hght is emitted, the slope of the composite flux power spectrum 
should therefore match the slope of the emissivity power spectrum. One might then expect a 
smooth roll-off from the slope of the emissivity power spectrum at these very low frequencies 
to a steeper slope at high frequencies, similar to what is observed in the hard states of GBHs 
and in AGN. In a simulation, the lowest reachable frequency is the inverse of the duration; 
unfortunately, the radius within which our simulation reached inflow equilibrium enclosed 
only a little more than half the luminosity, so we did not reach low enough frequencies to 
see the change in slope. 

One should also be aware of several other caveats in evaluating these comparisons with 
observations. First, dissipation in magnetically-dominated plasmas is thought to entail par- 
ticle acceleration across shocks or at reconnection sites. These processes can be much more 
rapid than the orbital timescale, while the inverse Compton cooling rate in units of f2x is 
~ {mp/me){Ls/ LE)~^{r / M)~^/'^ , for seed photon luminosity Lg. Thus, when there is signifi- 
cant thermal disk radiation, the inverse Compton cooling rate can likewise be much quicker 
than the orbital frequency. Because our cooling function has a characteristic rate ~ VLk-, it 
may underestimate high frequency variability. Second, our model focuses on the total coronal 
luminosity, which is likely dominated by photons at energies an order of magnitude higher 
than those generally studied in variability observ ations. If the power spe ctrum changes with 



photon energy (and there are some hints of this: iMarkowitz et al.l (120071 )). these may not be 
the appropriate comparisons to make. Third, our simulation did not distinguish the ther- 
modynamic properties of the disk body and the corona. It is possible that a more complete 
treatment of their thermal contrasts might alter the results. 

In conclusion, we have presented a radiative cooling model, based directly on simulations 
of 3D MHD turbulence in general relativity, that predicts the power spectra of fluctuations 
in hard X-ray flux observed from AGN and GBHs. The calculation used a new ray-tracing 
procedure for correctly tracking the propagation of light through time and space within the 
time-dependent 3D GRMHD simulation data set. The spectral slope found from our model — 
~ —2, significantly steeper than the slope of the accretion rate power spectrum, depends 
only weakly on the inclination and average accretion rate of the disk. Future simulations 
with more complete physics and a more complete traversal of parameter space will shed 
further light on this subject. 
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Fig. 1. — (a-d) Snapshots at — respectively — ip = {vr, 7r/2, 37r/2, 0}. Each image was taken at 
= 53° and t = 9000M with rh = 0.01. The axes mark the image's vertical and horizontal 
extent in the image plane in units of M. (Far Right) Logarithmic color map used to make 
the images. The intensity is normalized to the maximum intensity of the composite image 
shown in Figure O 




Fig. 2. — (Left) Composite of the snapshots shown in Figure [H (Right) Logarithmic color 
map used for the image. The maximum intensity in the map has value unity. 
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Fig. 3. — (Left) Normalized light curves (solid curves) for m = 0.01 and all values of 'd. The 
light curves and their mean values (dashed curves) have been shifted vertically by incremental 
factors of two for clarity. (Right) Normalized power spectra of these light curves compared 
to their best power-law fits (dashed curves). The power spectra are separated by incremental 
factors of 100. In both plots, the curves are ordered bottom-to-top in increasing order of 
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Fig. 4. — (Left) Normalized light curves (solid curves) for d = 29° and all values of m. The 
light curves and their mean values (dashed curves) have been shifted vertically by incremental 
factors of two for clarity. (Right) Normalized power spectra of these light curves compared 
to their best power-law fits (dashed curves). The power spectra are separated by incremental 
factors of 100. In both plots, the curves are ordered bottom-to-top in increasing order of rh. 
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Fig. 5. — Exponents a of the best power- law fits to the power spectra of the hght curves 
as functions of m (horizontal axis) and inclination angle (vertical axis). The departures 
from a ~ —2 in the upper-right-hand corner of the plot are caused by the disk's self- 
obscuration. The black curves represent contours of Ro- From bottom to top, Ro{^,m) = 
{risco,3.5M, 6M, 12M}. 
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Fig. 6. — (Left) Total flux F{t) observed at = 5° for m = 0.001 (solid black curve) 
compared to the accretion rate M (dotted blue curve) at r = 3.5M and emissivity C at 
r = 3.5M (dashed red curve). All rates are normalized to their time averages. (Right) 
Power spectra of these rates and their best-fit power-laws. The values of the best-fit power- 
law exponents are aj(_,j = —0.9, ac = —1.2, and ap = —2.1. 
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Fig. 7. — Ratio of the emissivity's normalized power spectrum to the accretion rate's normal- 
ized power spectrum plotted as a function of radius and frequency. Each power spectrum was 
smoothed over nine frequency bins before the ratio was taken in order to display trends in 
the data more clearly. Black curves show the orbital frequency (solid curve), radial epicyclic 
frequency (dashed curve), and vertical epicyclic frequency (dotted curve). 
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Fig. 8. — Ratio of the normalized power spectrum of dF/dr to the normahzed power spec- 
trum of the emissivity as a function of radius and frequency. Each power spectrum was 
smoothed over nine frequency bins before the ratio was taken in order to display trends in 
the data more clearly. 
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Fig. 9. — (Left) Phase ip^v^r) of dF/dr fluctuations when seen face-on (i? = 5°) for m = 
0.001. Note that we deviate from prior figure layouts and use a linear frequency scale here in 
order to resolve small-scale features. In addition, we show only the lower half of our frequency 
range. The dashed curve is the local inflow rate z/jnflow (Right) The linear, periodic color 
map used to generate this figure. 
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Fig. 10. — (Left) Normalized light curves obtained from the calculation ignoring time delay 
effects (dotted curve) and taking them into account (sohd curve). (Right) Normalized power 
spectra of these light curves compared to their best power-law fits; the dashed line represents 
the best fit to the data with time delay effects, the dash-dot to the data in which time delays 
were ignored. Both light curves are for = 29° and fa = 0.01. 
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Fig. 11. — Difference between the power-law exponents a from the calculation without 
time delays to that with time delays as a function of m (horizontal axis) and °d (ver- 
tical axis). The cases in the upper-righthand corner of the plot are heavily obscured. 
The black contour curves there represent — respectively — from bottom to top Ro{'&,'rh) = 
{risco,3.5M, 6M, 12M}. 



